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The interaction between natural selection and random mutation is frequently de- 
bated in recent years. Does similar dilemma also exist in the evolution of real net- 
works such as biological networks? In this paper, we try to discuss this issue by a 
simple model system, in which the topological structure of networks is repeatedly 
modified and selected in order to make them have better performance in dynami- 
cal processes. Interestingly, when the networks with optimal performance deviate 
from the steady state networks under pure mutations, we find the evolution behaves 
as a balance between mutation and selection. Furthermore, when the timescales 
of mutations and dynamical processes are comparable with each other, the steady 
state of evolution is mainly determined by mutation. On the opposite side, when 
the timescale of mutations is much longer than that of dynamical processes, se- 
lection dominates the evolution and the steady-state networks turn to have much 
improved performance and highly heterogeneous structures. Despite the simplicity 
of our model system, this finding could give useful indication to detect the underlying 
mechanisms that rein the evolution of real systems. 

I. INTRODUCTION 

The successful application of evolutionary process [l, 2] to the explanation of real networks, 
was one of the major achievements in the research of complex networks. Meanwhile, the 
topological structure of complex networks was frequently found to have determinative effect 
on the dynamical processes running on them^, |j, [5], and vice versajfj. Therefore, many 
dynamical systems were modeled as adaptive networks, in which 
processes was extensively coupled into the structural evolutions [ 



eedback from dynamical 
, y, 9]. In most of these 
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networks, the compositional elements could actively change their interaction subjects ac- 
cording to the dynamical states. But for many biological networks, structural mutation 
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y and randomly, with a timescale much separated from that of dynamical 
ll|. Besides, the structural evolution of these systems was mainly driven by 
preferential selection, in which systems exhibiting better performance during the life circle 
could pass their structural information to future generations with higher probability [l2|. 
Thus, an important issue emerged: how intensively could the pursuit of specific perfor- 



131 ]. reshape the topological 



mance, such as reaching functional states rapidly and robustly 
structures via preferential selection? 

In this paper, we investigate this question by designing a simple model system, in which we 
evaluate a network by its efficiency of escaping from disordered states in Local-Majority-Rule 
(LMR) dynamical processes(14| and try to evolve networks of low efficiency towards networks 



of high efficiency via repeated mutations and select ions [15|, [16|. Although LMR dynamics is 
too simple to represent most of the dynamical processes running on real systems, it could 
make the evolution only focus on optimizing the distribution pattern of edges among vertices, 
which greatly reduces the complexity of our problem. In the model system, we find the steady 
state of evolution depends strongly on the timescales of mutations and dynamical processes. 
When their timescales are comparable to each other, mutation dominates the evolution 
and the steady-state networks have similar structure as the steady-state networks of pure 
mutations. On the opposite hand, when the timescale of mutations is much longer than that 
of dynamical processes, selection dominate the evolution and highly heterogeneous networks 
with heavily connected hub and much improved efficiency emerge from the evolution. In the 
extreme situation, networks with optimal efficiency, which also deviate significantly from the 
steady-state networks of pure mutations in topological structure. At the end of this paper, 
we also propose a simple model on the evolution of particles. Additionally, we illustrate that 
this conclusion is still valid even in infinite population limit by a simple mathematical model. 
As an extension, this finding calls for a comprehensive understanding of the evolution of real 
systems, which is consistent with the suggestion that the effect of likelihood should be also 
incorp orated jl^] • 
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II. RESULTS 



In the evolution with low mutation rate (/i = 0.01), the networks evolve to have highly 
heterogeneous structures (Fig. dJ[a]), in which a vertex with degree comparable to network 
size N emerges to act as the communicating hub. Meanwhile, the degree distribution of the 
other vertices shifts to be power-law-like for small degrees, which is significantly different 
from the exponential degree distribution of the steady-state networks under pure muta- 
tion. Besides, these highly heterogeneous networks exhibit strong degree-degree correlations 
indicated by two parameters R ~ —0.495 ± 0.025, which means the global hub prefers inten- 
sively to interact with vertices having low degrees, and r « 0.067 ± 0.0083, which means the 
other vertices prefer to connect with vertices having similar degrees. We notice that such 
degree-degree correlations are consistent with the spin-spin correlations embedded in the 
strongly disordered spin configurations of highly heterogeneous networks. In other words, 
the dynamics-driven evolution is able to detect the correlation between different vertices' 
states and transfer it to topological clustering, which gives another potential explanation to 
the community-rich structures of many biological and social networks. 

On the opposite side, in the dynamics-driven evolution with high mutation rate (/i = 1) 
selection fails to achieve further improvement in both efficiency and structure than mutation- 
driven evolution, and the networks remain similar to the steady-state networks under pure 
mutations (Fig. HJ[b]). From extensive simulations, we find there exists a critical mutation 
rate fi c which classifies the dynamics-driven evolutions with different mutation rate into 
two distinct regimes (Fig. [2]). In the regime ji > /i c , mutations dominate the evolution, 
and the steady-state networks have similar structure as the steady-state networks under 
pure mutation. While in the regime \i < /i c , the steady-state networks appear to have 
highly heterogeneous structures, and their improvement in both efficiency and structure 
grows monotonically as fi decreased. Importantly, the critical mutation rate fi c decreases 
rapidly with N, which indicates a growing difficulty for large networks to achieve apparent 
optimization via dynamics-driven evolution. 

Interestingly, for the dynamics-driven evolutions with different sampling number Q, which 
determines the accuracy of efficiency sampling, similar transition emerges again (Fig. [3]): 



10 100 
Degree 



1000 



10 

Degree 



100 



FIG. 1: The steady-state network under dynamics-driven evolution with different mutation rate 
/i = 0.01 (a) or jj, = 1 (b) shows different vertex-degree distributions. Both of the networks are of 
size ./V = 2000, average connectivity (k) = 10 and minimal vertex-degree ko = 3. The solid lines 
represent the corresponding degree distribution for steady-steady networks under pure mutation 
(averaged over 200 samples), while the dashed line in (a) is the best power-law fit P(k) ~ i-92±o.02 
for the degree distribution of non- global hub vertices. 

only in the regime Q > Q c , highly heterogeneous networks emerge and their improvement 
grows monotonically with Q; besides, the critical sampling number Q c increases rapidly with 
N. Thus, we conclude that the dynamics-driven evolution is determined by the balance of 
two counteracting effects: the promotive effect generated from preferential selection strug- 
gles to push the population towards networks with high efficiency, while the degradative 
effect brought together by random mutation and imprecise sampling of efficiency drives 
the population back to the steady-state networks under pure mutations. Consequently, the 
steady state of dynamics-driven evolution should be independent of the original networks, 
even when they are optimally designed (Fig. H]). 

To give a mathematical illustration, we consider a population of M particles with two 
possible energy states S\ = and 62 = 1. For a particle with energy £j (i = 1 or 2), it will 
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FIG. 2: Transition of the steady-state mean efficiency £ , mean hub-spoke coefficient Ch s , mean 
correlation index R and mean assort ative-mixing index (of global hub-removed subnetwork) r under 
dynamics-driven evolution versus mutation rate ft. The critical mutation rate [i c = 0.25, 0.2, 0.15 
for network size N = 500 (diamonds), 1000 (circles), 2000 (squares), respectively. Each point 
represents a average over 6 x 10 4 generations after the steady state of evolution is reached. 



be assigned with an random evaluation following normal distribution 



(1) 



Q <Ji V 27T 

here the standard deviation coefficient Oi satisfies 1 <S of <g; +oo and Q e [l,+oo) is the 
control parameter. Besides, the offspring particle generated from a particle with energy 
will appear on the opposite energy level with probability Pi, otherwise it will have identical 
energy with its parent. For mutation- driven evolution, every particle is replaced by an 
offspring generated from it at each generation. Easily seen, the particle composition r g , 
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FIG. 3: Transition of the steady-state mean efficiency 8, mean hub-spoke coefficient Chs, mean 
correlation index R and mean assort ative-mixing index (of global hub-removed subnetwork) r under 
dynamics-driven evolution versus sampling number f2. The critical sampling number Q c = 220, 
300, 420 for network size N = 500 (diamonds), 1000 (circles), 2000 (squares), respectively. Each 
point represents a average over 6 x 10 4 generations after the steady state of evolution is reached. 



which is defined to be the fraction of particles with energy e 2 at generation g, will finally 
evolve to a steady state 

Pi 



(2) 

Pi +P2 

While in evaluation-driven evolution, every particle generates E offsprings (E is positive and 
finite) at the start of each generation, then only the M particles with the highest evaluation 
of all (E + 1)M particles will survive and pass into next generation. Therefore, in the infinite 




FIG. 4: (Color Online) The evolution of mean efficiency £, mean hub-spoke coefficient C^s, mean 
correlation index R and mean assortative-mixing index (of global hub-removed subnetwork) r as 
a function of simulation generations. The mutation rate is fi = 0.05 (black and green curves) or 
fi = 0.01 (red and blue curves). In the dynamics-driven evolution started from hub-spoke networks 
(black and red curves), the population is refreshed mostly by degraded networks due to their 
dominant fraction in the offsprings, and gradually converged the steady state of dynamics-driven 
evolution started from random Poissonian networks (green and blue curves). 

population limit (M — > oo), the steady-state particle composition r| will satisfy 

r e = [r e + E{1 _ r e )pi + ^ _ ^ f+°°^^ °A) dx 

1 - rt = [1 - rt + E(l - rt)(l - Vx) + Er e s p 2 ] Mfa, §)dx ' 

here <p s is the steady-state threshold that only particles with evaluation higher than it could 
survive. When Q — > +oo, which means the stochastic fluctuation in particle evaluation 
will decrease to vanish, we find that J^°° N{e\, ^)dx — > J^°° N{s\, 0)dx = and r| — > 1. 
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Oppositely, when f2 is sufficient small so that 

f+oo 2 r+co 2 

/ tf( £l , jr)dx « / tf(e a , ?*)ds, (4) 

which implies particles with different energy will survive with approximately equal proba- 
bility, the steady-state composition will have r% ~ Pl + P2 = r T- O n the other hand, in the 
extreme condition pi — > and p 2 —>■ I, which corresponds to the dynamics-driven evolution 
with high mutation rate, it could be derived from equation (j2J) and equation (j3J) that r e s — > 
and r™ — > 0. In other words, the steady-state population of evaluation-driven evolutions has 
similar composition with that of mutation-driven evolution. However, in the extreme condi- 
tion pi ?a p 2 <C 1, which corresponds to the dynamics-driven evolution with small mutation 
rate, we find that r\ ~ 1 and r™ ~ ^. Obviously, this compositional difference indicates 
an apparent optimization achieved by evaluation-driven evolution. Thus, if we take the 
dynamics-driven evolution as particles jumping among a series of structural heterogeneity 
levels, it could be easily derived from above analysis that the optimal structure could be 
approached only when the timescale of mutations is much longer than that of dynamical 
pro cesses jl8|]. 

Finally we replace the local rewiring scheme by preferential rewiring, and show there also 
exists a critical strength /3 C (Fig. [5]) to classify dynamics-driven evolutions with different 
strength of preferential rewiring into mutation- dominating regime (/3 < /3 C ) or selection- 
dominating regime (/3 > /3 C ). Interestingly, the critical strength f3 c is found to increase 
monotonically with iV as well, which reemphasizes the growing difficulty to achieve apparent 
optimization in large networks and supports that our finding is independent of the way how 
mutation occurs. 



III. DISCUSSION 



It has been debated whether the evolution in biological world should be viewed as some 
hill-climbing processj^], or is mainly determined by the "neutral mutations" [21I] like a 
random walk. But in our simple model system, we find when the optimal networks devi- 
ate from the steady state networks under mutations, a flexible combination of these two 
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FIG. 5: Transition of the steady-state mean efficiency £ , mean hub-spoke coefficient Chs, mean 
correlation index R and mean assort ative-mixing index (of global hub-removed subnetwork) r under 
dynamics-driven evolution versus the strength of preferential rewiring f3. The critical strength /3 C = 
0.9, 0.95, 0.99 for network size N = 500 (diamonds), 1000 (circles), 2000 (squares), respectively. 
Each point represents a average over 6 x 10 4 generations after the steady state of evolution is 
reached. 

viewpoints is required. Interestingly, similar deviations are frequently detected in biologi- 
cal and social systems, such as the abundance of certain network motifs 



22 or nucleotide 



types[23]. Therefore, dynamics-driven evolution could be a potential way to explain these 
structural heterogeneities. However, the evolution of real systems is far more complicated 
than our simple model system. Meanwhile, the evaluation of them is usually determine d by 
a combination of multiple components such as efficiency, sensitivity, robustness and cost [24 ] . 
Therefore, more detailed models are still needed to give a comprehensive understanding to 
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the evolution running in real world. 



IV. METHODS 



A. Local-Majority-Rule dynamics 



For a undirected network of size N and average connectivity c, we attach to each vertex i 
a binary spin variable G{ — ±1. At each time step t of the LMR dynamics, all the vertices 
update their states synchronously according to 

«T i (t)=sigIl[^(T i (t-l)], (5) 

jedi 

here di means the vertices which are connected with vertex i directly. As described 

n 

in [181 ]. we start a LMR dynamical process from a strongly disordered configuration 

N 

a(0) = {o"i(0), (72(0), ... , ctat(O)}, which is randomly generated with the constraints ^ Oi = 

i=l 

N 

^2 ki<Ji = (ki is the degree of vertex 1). Typically, the LMR dynamics will drive the network 
i=i 

to stabilize at a consensus state, in which all the vertices have the same spin state. It has 
been shown in [l^] that in LMR dynamical processes, scale-free networks with degree distri- 
bution P(k) ~ fc -7 and 7 — > 2 can reach consensus more rapidly than Poissonian networks 
of similar size and average connectivity. In other words, highly heterogeneous networks are 
more efficient in communicating internal states than homogenous networks. 



B. Definition of efficiency 



For a given network Q, we start a total number of Q different LMR dynamical processes, and 
run them 
shown 

is determined by the the escaping velocity of the network's spin configuration from the 
strongly disordered region. In this paper, we define the efficiency Q of network Q as the 



"or T = 1 time step to reach the corresponding configuration a a (l). It has been 
in lj| that the characteristic relaxation time of a network in the LMR dynamics 
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average change of vertex-state in <?(1), which is calculated according to 



U N 



(6) 



a=l i=l 



Easily seen, the efficiency sampled in this way is only a random estimation of the true 
efficiency E*(Q) = lim S(Q). According to the Central Limit Theorem, the distribution 



of £{G) will be approximately Gaussian-distribution-like with the average value equal to 
£*{Q)i while the mean square variance will be proportional to 1/Q (see the Supplementary 
Information). That's to say, the efficiency of a network will be evaluated more precisely with 
large Q. 



The dynamics-driven evolution starts from a population of P networks, which are uniformly 
sampled from the ensemble of random Poissonian networks of given size N and average 
connectivity (k). At the beginning of each generation, every network generates E exact 
copies, which expands the population to a size (E + 1)P. Then every new network undergoes 
random mutations following local rewiring [251]: 




(1) For each vertex, mutation occurs with probability fi (mutation rate); 

(2) if mutation is accepted by vertex i, a rewiring of edge — > is proposed, in 
which vertex j is randomly selected from the nearest neighbors of vertex i and vertex 
I is then randomly selected from the nearest neighbors of vertex j except vertex i; 

(3) this proposal will be rejected if and only if edge (i, I) already exists or the degree of 
vertex j is less than a minimal value ko after cutting edge 

Then the efficiency of each network is sampled independently, and the population shrink to 
its original size with only P networks with the highest efficiency survive and pass into next 
generation. As a reference to the dynamics-driven evolution, we also derive the steady-state 
networks under pure mutations by proposing a mutation- driven evolution. In the evolution, 



C. Evolutionary algorithm 
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all the networks keep to be replaced by a new network generated from it following the same 
rules of mutation. For the local rewiring scheme, the degree distribution of the steady-state 
networks under pure mutations is shown to decay exponentially 3] . If not specified by the 
context, the parameters default to adopt value Q = 1000, P = 25, N = 1000, (k) = 10 , 
E = 3 and k = 5. 



In this work we also involved preferential rewiring scheme [26|]: the only change from above 
description is that vertex I (I ^ i and / ^ j) in the rewiring proposal is chosen randomly 
with a probability Hi oc (ki + l) 13 , here k\ is vertex Z's degree and (3 is the control param- 
eter which determines the strength of preferential rewiring. Obviously when (5 = 0, the 
preferential rewiring scheme will degrade to random rewiring, under which the steady-state 



networks has been proved to be Poissonian 2g] . On the opposite side, the maximal value 



of P used in this paper is 1, under which the steady-state networks will have exponential 
degree distributions [ 



D. Identification of global hub and Definition of hub-spoke coefficient 

For generality, we identify the vertex with the highest degree in a network as the global hub 
for both heterogeneous and homogenous networks. (If there are multiple vertices with the 
highest degree, we only identify the one with the smallest index as the only global hub.) We 
also define the hub-spoke coefficient Chs of a network to be Chs = lr-i> i n which khub is the 
global hub's degree and N is network size. Obviously, when Ch s equals to 1, the network 
will look completely hub-and-spoke like. 



E. Quantification of degree-degree correlation in networks 

To give an explicit quantification of the degree-degree correlations in a network, we define two 
quantities: correlation index R measures the degree-degree correlation on links connecting 
the global hub and the other vertices, while assortative-mixing index r measures the degree- 
degree correlations on links connecting the other vertices. We define R to be the normalized 
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R 



ratio of the mean degree of nearest neighbors of the global hub (k^) to the averaged 
value (knn b )ran of this mean degree over an ensemble of randomly shuffled networks [221] . and 
calculate it according to 

(/UhubX _ jh.hub\ \ I ( / UhubX _ / UhubX \ :f I h.hub\ ^ / L,hubX 
XX^nn I X^nn Iran)! KX^nn I max X^nn Iran)) 11 X^nn I ^ X^nn Iran 

(lh,hubX _ /L,hub\ \/(/t,hub\ . _ / UhubX \ ■£ / UhubX ^ / UhubX (J) 
XX^nn I \™nn Iran) I \\^ nn I mm X^nn Iran)-, 11 X^nn I X^nn Iran 1 \' I 

n ;f lh,hub\ _ I h,hub\ 
u i 11 X^nn I X^nn Iran 

here (k!^) max and (&^ 6 ) min are the maximal and minimal value of (fc^f ) in the ensemble 
of randomly shuffled networks, respectively. Obviously, R G [—1,1]. When R < 0, the 
global hub prefers to interact with vertices with low degrees compared with its behavior in 
the randomly shuffled networks. Otherwise, when R > 0, vertices with high degrees are 
preferred. On the other hand, we define r to be the assortative-mixing index of the globe 
hub-removed subnetwork following [27| 

' M' 



M' 
i=l 



E 

i=l 



jj-\-kj 
2 



M' 



i=l 



M' 

E 

i=l 



2 



2 • 



(8) 



here M' is the number of edges in the subnetwork. When r > 0, these vertices prefer to 
be assortatively connected with vertices having similar degrees; Otherwise, when r < 0, the 
subnetwork is disassortatively connected. 



To make rapid calculation of R, here we give a theoretical prediction of ) mn by 



lh hub \ 
\ K nn I 



k 



— ^k(N k -5 k , kh JP*(k), 

hub k=1 



(9) 



in which N k is the number of vertices with degree k, 5 is the Kronecker symbol and P*(k) is 
the probability that a vertex with degree k is connected with the global hub in a randomly 
shuffled network. Now we approximate P*(k) by 

p*(k) « £ \ p ^ k ) ID 1 - > ( 10 ) 



8=1 



3=1 



in which 



2M - k hub -EE W(N k , - 5 k ., kh JP*{k>) EI (1 - P m (k'))} 

j=l ^fc'=l m=l 
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In numerical simulations, it's shown to be a good approximation for the networks appeared 
in this paper (see the Supplementary Information). 



F. Generation of hub-spoke networks 



A hub-spoke network of size N, average connectivity c and minimal degree ko is generated 
through: 

(1) first generate a random Poissonian network of size N — 1 and average connectivity 

jV-1. 



c = c 



N ' 



(2) then evolve the Poissonian network under mutation-driven evolution with minimal 
degree k' Q = k — 1 until the steady state is reached; 



(3) then assortatively shuffle the mutated network following 28] with p = 1 



(4) finally add a global hub to the shuffled network and connected it with all the other 
N — 1 vertices. 

From numerical simulations, we find such hub-spoke networks have much higher efficiency 
than the Poisson networks of the same size and average connectivity (Fig. H]). 
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I. LOCAL-MAJORITY-RULE DYNAMICS 

A strongly disordered configuration could be constructed as follows: one first divide the 
vertices into two groups (G+ and GJ) of equal size randomly and assign spin +1 to vertices 

N 

in group G + and assign spin —1 to vertices in group G_; then as long as ^ fc^cr, ^ 0, two 

i=i 

vertices (one from G + and the other from GJ) are chosen randomly and their states are 

N 

exchanged if and only if this exchange does not cause an increase in the value of | ki a i I ■ 

i=i 

In such a way, a vertex (either chosen randomly or reached by following a randomly chosen 
edge) will be found to be in the spin +1 state with probability 1/2. 

Easily seen, there exist clear spin-spin correlations in the strongly disordered configurations 
of highly heterogeneous networks. For a network, we define the probability distribution q(k) 
as the probability that a vertex with degree k is found to have the same spin state as the 
global hub in the strongly disordered configurations. Then we find for highly heterogenous 
networks, q(k) is clearly higher than 1/2 when k is small and is clearly lower than 1/2 when 
k is large (Fig. ??). That's to say, the degree-degree correlations in the highly heterogenous 
networks reached by dynamics driven evolution are consistent with the spin-spin correlations 
embedded in the strongly disordered configurations of these networks. 

However, these spin-spin correlation aren't the main causation of the emergence of highly 
heterogenous networks in the dynamics-driven evolution. To prove it, we modify the model 
system by using simplified LMR dynamics. In the simplified LMR dynamical processes, 
the initial spin configurations are disordered configurations which only need to satisfy the 

N 

constraint ^2 o% — 0. In this paper, a disordered configuration is constructed by dividing 
i=i 

the vertices into two groups of equal size randomly and assigning opposite spin states to 
the vertices in different groups. From numerical simulations, we find the dynamics-driven 



